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Abstract 

Within the theory of Siegert pseudostates, it is possibie to accurately calculate bound states 
and resonances. The energy continuum is replaced by a discrete set of states. Many questions 
of interest in scattering theory can be addressed within the framework of this formalism, thereby 
avoiding the need to treat the energy continuum. For practical calculations it is important to 
know whether a certain subset of Siegert pseudostates comprises a basis. This is a nontrivial issue, 
because of the unusual orthogonality and overcompleteness properties of Siegert pseudostates. 
Using analytical and numerical arguments, it is shown that the subset of bound states and outgoing 
Siegert pseudostates forms a basis. Time evolution in the context of Siegert pseudostates is also 
investigated. From the Mittag-Leffler expansion of the outgoing-wave Green's function, the time- 
dependent expansion of a wave packet in terms of Siegert pseudostates is derived. In this expression, 
all Siegert pseudostates — bound, antibound, outgoing, and incoming — are employed. Each of these 
evolves in time in a nonexponential fashion. Numerical tests underline the accuracy of the method. 

PACS numbers: 03.65.Nk, 02.70.-c 
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I. INTRODUCTION 



Let us consider the following model potential: 

f -V , < r < a 
V(r) = { ° ' - , (1) 

I , r > a 

where Vq > and a > 0. In order to find real-energy solutions, one normally matches 
the wave function inside the well to a superposition of incoming and outgoing plane waves 
outside the well (or to an exponentially damped function in the case of a bound state). 
Following Siegert Q], however, we require that the solution for r > a is proportional to 
exp (ikr), i.e., we apply the Siegert boundary condition 

d , 



dr 



\kip{a) . (2) 



If, in addition, we demand vanishing of the wave function at the origin, the wave number k 
must satisfy the following relation: 

ik = \/k 2 + 2V cot (v 7 ^ 2 + 2V Q a) . (3) 

This transcendental equation can be solved for only a discrete set of k values. As was 
done, for instance, in Ref. Q], we solve Eq. @ numerically, after setting Vo = 5 and 
a = 10. The resulting discrete k spectrum is presented in Fig. There are 10 bound 
states, which appear on the positive imaginary axis. The nine antibound states lie on the 
negative imaginary axis. The solutions with Re (A;) > are associated with outgoing Siegert 
pseudostates. Similarly, states with Re(fc) < refer to incoming Siegert pseudostates. 

Traditionally, Siegert-state theory has focused on scattering resonances and decaying 
states. See Refs. BQB0 for an overview of the literature. Several methods exist that 
allow one to directly calculate the complex energy of a decaying state: complex scaling 
DQD , complex absorbing potentials [lOj, |ll| , and the direct solution of the Schrodinger 
equation subject to the Siegert boundary condition, Eq. (J2J). This third approach used to 
be numerically inefficient, because a nonlinear eigenvalue problem had to be solved in an 
iterative fashion LiLiLi- 

nn 

A major step forward was made by Tolstikhin, Ostrovsky, and Nakamura p, Their 
method provides, after solving a single generalized eigenvalue problem, access not only to 
the bound, antibound, and resonance states, but also to the discretized pseudocontinuum. 
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FIG. 1: The complex k spectrum for a particle in the step potential of Eq. Q with Vq = 5 and 
a = 10. The spectrum is entirely discrete. Ten bound states [Re(fc) = 0, Im(fc) > 0] are present, as 
well as nine antibound [Re(fc) = 0, Im(fc) < 0]. The positive Re(fc) branch of the spectrum shows 
the complex wave numbers of the outgoing Siegert pseudostates, the negative Re(fc) branch the 
incoming. 

In addition, it becomes possible to derive fundamental properties of Siegert pseudostates 
using simple and elegant mathematical techniques. An application of the method of Refs. 

|J to a molecular fragmentation problem is the subject of Ref. [3]. In that work, only 
bound states and outgoing Siegert pseudostates were utilized as fragmentation-channel basis 
functions. 

A point of central concern in this paper is the time evolution of a wave packet expanded 
in terms of Siegert pseudostates. This aspect was first addressed by Yoshida et al. [16]. 
They successfully introduced Siegert pseudostates as a basis capable of eliminating artificial 
boundary reflections. In a subsequent paper [17], it is argued that the time evolution for 
t > is given by 



where the imaginary part of the complex energy E n is negative. f n { r ) is the spatial repre- 
sentation of the nth Siegert pseudostate, and 




(4) 



n 




(5) 



In Sec. |nj we review some ideas of the formalism of Tolstikhin, Ostrovsky, and Nakamura. 
We extend their work by providing arguments why certain subsets of Siegert pseudostates 
may be employed as bases. Moreover, we investigate time evolution in the context of Siegert 
pseudostates. We find that rigorous Siegert-pseudostate theory, in general, necessitates 
nonexponential time evolution, in contradiction to Eq. (J3J). Numerical evidence is presented 
in Sec. IIHI The calculations demonstrate that for a fast-moving wave packet Eq. (£Q) is 
accurate, while its performance deteriorates as the average energy of the wave packet is 
decreased. Section |IV] concludes. Atomic units are used throughout. 



II. MATHEMATICAL CONSIDERATIONS 

Consider the radial Schrodinger equation 

H<p(r) = E V {r) , (6) 

where 

V(r) is an arbitrary (short-range) potential. We seek solutions </?(r) for r e [0, a] (a > 0), 
such that at r = the usual boundary condition for the radial eigenvalue problem 

<P(0) = (8) 

and at r = a the Siegert boundary condition, Eq. (J2J), are satisfied. The solutions <p(r) are 
referred to as Siegert pseudostates jsHsj]. Equation (J2J) allows us to separate solutions to the 
Schrodinger equation into purely outgoing, purely incoming, bound, and antibound states. 

The wave number k in Eq. (J2J) is in general complex. Its relation to the eigenenergy E is 
given by 

E = y + V(a) , (9) 

which implies that a should be chosen sufficiently large, if possible, such that V(r) = const, 
for r > a. Applying Eq. Q in the case of a long-range potential introduces the approxima- 
tion of forcing the potential to be constant beyond r = a, while preserving continuity of the 
potential at that point — in contrast to the approach taken in Ref. P, 0], which enforces a 
discontinuity by setting V(a) = in Eq. (jHJ). 
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Utilizing a set of N linearly independent, not necessarily orthogonal basis functions 
{yj(r) : j = 1, . . . , N} on the interval [0, a] and assuming completeness in the limit iV — > oo, 
we can make the ansatz 

N 

<f(r) = ^2cjyj(r) ■ (10) 

i=i 

If we insert this into Eq. ©, multiply from the left by yi{r), and integrate over r from to 
a, we find 



~o I ^( r )T^yi c i%( r ) dr + / yi( r ) V ( r )y2 c jyj( r )dr = E / ^(r) V c j ?/ i (r)dr (11) 
2 Jo &r . =1 J j=i J . =1 

The term containing the second derivative with respect to r can be integrated by parts. 
Thus, after applying the boundary conditions in Eqs. © and (0), we define 



H *i = / d^( r )cV^ r ) dr (12) 



o 



+2 / yi{r)[V(r)-V(a)] yj (r)dr 



s a = / yi( r )yj( r )dr , (13) 

Jo 

Lij = yi(a)yj(a) , (14) 

h — ik , (15) 

and arrive at the nonlinear eigenvalue problem 

H + h 2 S - k£) c = . (16) 



We refer to the eigenvector c in this equation as Siegert pseudovector. In practical calcula- 
tions, the basis functions yj(r) may be chosen to be real. Thus, the matrices H, S, and L 
may be assumed to be symmetric and elements of M 7Vx7V . Note, however, that the vector c, 
consisting of the expansion coefficients Cj, j = 1, . . . , N [see Eq. (J10J) ]. is complex in general: 
c G C . If Eq. f!16|) is satisfied by an eigenpair (x, c), then — in view of the realness of H, 



S, and L — the pair fx* c*) also forms a soi 

nn 

Tolstikhin, Ostrovsky, and Nakamura p, |6| 



ution. 

observed that the nonlinear eigenvalue prob- 



lem of Eq. (jlfij) can be recast in terms of a linear one by defining 

d = xc . (17) 
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Then, Eq. (ITtjj) goes over into 



-Hc= x(Sd-Lc) 



and, trivially, 



Sd = xSc . 

These two matrix equations are equivalent to the generalized eigenvalue problem 

Ax = kBx . 

The symmetric matrices A and B 6 ~^Nx2N are gj ven \yy 

H 

A 

S 



(18) 



(19) 



(20) 



B 



-L S 
S 



(21) 



(22) 



and the eigenvector x e C 2N in Eq. (J2*U|) is composed of the two vectors c and d: 



x 



c 
d 



(23) 



Let us denote by N the number of linearly independent solutions x n (eigenvalue x n ) to 
the generalized eigenvalue problem in Eq. ()20|). (It is by no means clear that iV = 2N.) We 
then define the matrices 

I = K...,^]eC 2 ™ (24) 

and 

K = diag(x 1 ,...,x J v)GC^, (25) 

such that 

AX = BXK . (26) 

The overlap matrix S [Eq. 1)13)1 ] is positive-definite and therefore invertible. Its inverse, 
is also symmetric. It can be concluded that the inverse of B exists, since it is easily 
seen by direct construction that 



B 



-i 



5 1 

S- 1 S- l LS- 1 



(27) 
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The real symmetric matrix B can be diagonalized via an orthogonal similarity transforma- 
tion, 

U T BU = D , (28) 

where U E R 2Nx2N i s orthogonal, U T U = 1, and D E R 2Nx2N i s diagonal. All diagonal 

elements of D differ from zero, for D is similar to B and thus also invertible. With this in 
mind, the matrix 

U = UD- 1 ' 2 E C 2Nx2N (29) 

can be introduced in a meaningful way. While B is invertible [Eq. ([27))]. it is in general 
neither positive- nor negative-definite. Hence, some of the columns of U are real, but the 
others are purely imaginary. The matrix U can be utilized to convert the real symmetric, 
indefinite generalized eigenvalue problem, Eq. (|2T)j) . to a standard one: 

AX = XK . (30) 

Here, 

A = U T AU (31) 

and 

X = DU T X . (32) 

The matrix A is complex symmetric. Unfortunately, complex symmetry is not a very 
useful property, as it is known that any complex matrix of square format is similar to a 
complex symmetric matrix (see, for instance, Ref. [10]). Generally, it is not possible to 
guarantee diagonalizability, i.e. the existence of a basis of eigenvectors, of a complex sym- 
metric matrix. Note that A consists of purely real and purely imaginary matrix blocks, but 
whether this helps to prove its diagonalizability is currently unclear. A sufficient condition 
for diagonalizability is that all 2N eigenvalues x n are distinct. Degeneracies could cause 
difficulties, but in numerical calculations true (accidental) degeneracies are practically never 
encountered. Under the assumption that A is diagonalizable, N = 2N and the matrix 

X=[x 1 ,...,x 2N ]EC 2N * 2N (33) 

may be chosen to be complex orthogonal, 

X T X = 1 , (34) 
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as shown, e.g., in Ref. 10]. The 2N vectors x n are linearly independent and, consequently, 
form a basis of C 2N . This fact can be expressed in compact form in terms of the completeness 
relation 

2N 

$>n5£ = l. (35) 



71=1 



From Eqs. an d it follows for the solutions x n of the original generalized eigen- 
value problem in Eq. (|20|) that 



x^Bxn = S mn , m, n = 1, . . . , 2N , 



and 



2 at 

71=1 



-1 



(36) 



(37) 



Combining Eq. (j3*oT) with Eqs. (|17|h (|22j). and (|23j). and employing the normalization 
convention of Refs. [l], Q|, the orthonormality relation satisfied by the eigenvectors c n of the 
nonlinear eigenvalue problem in Eq. (J16j) reads 



'-m k -"-n 



ciLc n 



&mn i TTlj TL 1, 



2N 



(38) 



Take notice of the occurrence of a singularity in this expression if x m + x n approaches for 
some pair of indices m, n. This happens whenever there is a deeply bound state (eigenvalue 
x m ), since there exists a corresponding antibound state whose eigenvalue x n equals — x m to 
machine precision. 

Even in the absence of the term involving the surface matrix L, Eq. (J38J) could not be 
used to establish a proper inner product on C^, because, for a given vector v e C N \ {0}, 
v T Sv is not necessarily different from zero. The occurrence of indefinite inner products 
is also well known in the context of complex scaling [3, 0] and in the method of complex 
absorbing potentials jll], 
fl23J), (E7J, and (JUT}: 



111 ]. The following three relations are consequences of Eqs. (fTTjl . 



V —c n cl = , 
z — ' x„ 

71=1 " 
2N 

C n C n — 1 , 

n=l 

^x„c„c^ = 2S- 1 LS- 1 . 



(39) 
(40) 
(41) 



n=l 



The normalization suggested by Eq. pa}) has been applied. 

Equation (|3Uj) demonstrates that any vector v G C N can be represented as a superposition 
of the Siegert pseudovectors c n : 

1 2iV 

« = 2 E ( c n Sv ) c « • ( 42 ) 

n=l 

This representation, however, is not unique, for the vectors c n , n = 1, . . . , 2 A, form an 
overcomplete subset of C^. The rank of the matrix C — [c 1; . . . , c 2 n\ eC Nx2N cannot be 
greater than A [in fact, because of Eq. (jjUj) , rank(C) = A]. Thus, the vectors c 1; . . . , c 2 n 
are linearly dependent, i.e., the equation 

2JV 

^« n c n = (43) 

n=l 

can be satisfied by some a n 7^ 0. If we define 

M mn = clSc n , (44) 

then 

2N 

J2M mn a n = ,m = l,...,2N , (45) 

n=l 

from which we may conclude that the matrix M G C 2Arx2JV cannot be invertible (since not 
all a n have to be equal to 0). The nonuniqueness of the representation in Eq. (J4*2|) is directly 
linked to the noninvertibility of M. 

Another instructive way of looking at this is to make use of Eq. (|4f)|) to derive the 
following relation: 

MM = 2M . (46) 

If M were invertible, then we would have M = 21. This contradicts the orthonormality 
relation, Eq. (jSSjl: 

c T Lc 

]\/f — X I ^rn-^^n ( ±7} 

%m ~r l^n 

For example, for m = n the surface term can be written as ip m (a) 2 /2x m [Eqs. (110)1 and 
(|14|) ]. which differs from unity in general. 

In Appendix |X} we demonstrate that a subset of {ci, . . . , C2at} can be found, comprising 
exactly vectors and forming a basis of C^. Without loss of generality, the elements of 
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this subset are taken to be the vectors Ci, . . . , Cat. The proper completeness relation allows 
one to represent an arbitrary vector v G C N in a unique way: 

N 

v = E a n c n . (48) 

n=l 

Necessarily, 

N 

J2 M m n a n = clSv ,m = l,...,N . (49) 

n=l 

Note that the matrix elements M mn in this equation refer only to the selected subset. The 
linear independence of the Siegert pseudovectors c 1; . . . , c N ensures the invertibility of M G 
C NxN [cf. Eqs. (jUIJ), (HU), and (jUJ)]. The expansion coefficients a m are unique and are 
given by 

a m = ^(M- 1 ) mn c^. (50) 

n=l 

The associated completeness relation reads 

Af Af 
m=l n=l 

Let us now turn our attention to the problem of describing the time evolution of an 
outgoing wave packet utilizing Siegert pseudostates. For that purpose, it is natural to 
choose for the basis {ci,...,Cjv} all outgoing Siegert pseudovectors [Re(/c n ) > 0] plus a 
complementary number of bound eigenvectors. (We provide in Sec. IIHI numerical evidence 
that these form indeed a basis of C .) Consider an initial wave packet ip(r;t = 0) that is 
entirely localized within the interval [0, a]. Using Eqs. ()10|) and (|51|) . and assuming pure 
exponential time evolution for each Siegert pseudostate, the time evolution of the wave 
packet for t > is obtained from 

Af Af 

^ *) = E E ( M_1 Ln WVOe^Wr) , (52) 

m=l n=l 

where 

£ m = ^ + V(a) . (53) 

One disadvantage of Eq. (|52j) is the need to invert a matrix. In addition, as we will see below, 
the assumption of exponential time evolution is wrong in general for Siegert pseudostates. 

A practically and formally more acceptable expression for the time evolution of the wave 
packet can be derived from the Mittag-Leffler partial fraction decomposition 3| of the 
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outgoing-wave Green's function represented with respect to Siegert pseudostates. This rep- 
resentation has been given in Refs. J, 0,0]. The outgoing- wave Green's function satisfies 
the equations 



(E — H)G(r, r'; k) = 5{r-r') , 
(7(0,7"'; k) = , 
—G(r,r';k) = ikG(a,r';k) 

for r, r' G [0, a]. These conditions can be translated to 



(54) 
(55) 
(56) 



H + x 2 5 - xL) G{k) = -21 , 



(57) 

where the matrix representation of the outgoing-wave Green's function in the basis of the 
functions Vj{r) is defined via the relation 

N N 



G(r,r'; k) = £ £ G tJ (k) yi (r)y 3 (r') . 

i=l j=l 



Making the ansatz 



2iV 



G(k) = ^ a nC n cl 



(58) 



(59) 



n=l 



and putting Eqs. (|39|). (|40p to use, it follows from Eq. (J57)l that 

1 



Hence, 



a, 



2iV 



k n (k k n ) 



G(r,r';k) = J2 



<p n (r)<p n (r / ) 



r, r' G [0, a] 



(60) 



(61) 



A; n (A; /c n ) 

n=l v ' 

serves as an approximation to the outgoing-wave Green's function, within the framework of 
the underlying finite basis set {yj(r) : j — 1, . . . , iV}. 

The Green's function in Eq. (pTTj) allows us to determine the time evolution of the wave 
packet for t > (l9|: 



l 



2N 



Hr;t) = — I dEe~ m dr'G(r,r'; k)^(r';t = 0) = V/M*)^)^ 



(62) 
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Here, 



2tt 



-dE 



(63) 



The integration path lies on the physical sheet [3J of the energy Riemann surface, and runs 
from — oo to +00 infmitesimally above the real E axis (lflj . In order to evaluate /3 n (t), we use 
Eq. and replace the integration path in the energy plane by the integration contour in 
the complex k plane illustrated in Fig. Et- For t < 0, this contour can be closed in the first 
quadrant of the k plane without affecting (3 n {t). The resulting integration loop is indicated in 
Fig. 12b. The integration path along the positive imaginary axis is infmitesimally displaced 
toward positive real k values, so that the integrand in Eq. ()63)1 is analytic on and inside the 
loop. Hence, (3 n {t) = for t < 0, as expected for a description based on the outgoing-wave 
Green's function. 

If t > 0, then (3 n (t) can be written as the sum of two separate contour integrals. The 
respective contours are shown in Fig. Et- Thus we have 



(64) 



where 



and 



e lEnt for bound states 
otherwise 



i exp {— iV(a)t} 
2n k n 

exp {— iV(a)t} 

g-iTT/4 fa 



k 



k kn 



-e- ikt / 2 dk 



2tt 



kn 




t 



00 e -ifc 2 t/2 
k — h 



dk 



From Eq. (|39j) it follows that 



2N ^ 



fr) = 



n=l 



Therefore, the term proportional to t in /% (t) does not contribute to r/>(r; t), Eq. (jB^j) 



(65) 



(66) 



(67) 



and we will not consider it further. As demonstratec. 
can be expressed utilizing the Faddeeva function 



in Appendix iBl the integral in Eq. 



20 



21 



22, 



231,124, 



25 



w(z) 



-ds , Im(^) > 



(68) 
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a) 



lm(k) 



Re(k) 




c) 



lm(k) 



Re(k) 



lm(k) 



Re(k) 



FIG. 2: a) The integration contour in the complex k plane employed to evaluate Eq. (jfiMj) . The 
integration proceeds along the imaginary axis from ioo to the origin, and from there along the real 
axis to +00. This contour corresponds to integration from —00 to +00 on the physical sheet of the 
energy Riemann surface, b) If t < 0, then the contour may be closed in the first quadrant of the k 
plane, c) The integration path in a) may be evaluated for t > using one closed contour and one 
that runs from —00 to +00. 
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Hence, on combining Eqs. (JMJ), (J55j) . and the results from Appendix iBl the time- 

evolution coefficients /3 n (t), t > 0, read 



e lfc ™*/ 2 — |iy fe 17r / 4 y|A; n ^ for bound or outgoing states 



/?„(t) = exp{-iU(a)t} x < 



~w e m ^\J\k r ^j for antibound or incoming states 

(69) 



(The square root in the argument of the Faddeeva function must be interpreted as a positive 
number.) The result in Eq. (|69|). together with Eq. (|62jh represents — within the framework 
of Siegert pseudostates — the most consistent approach to describing the time evolution of 
a wave packet. It is interesting to observe that in the limit t — > + , all j3 n (t) tend to 1/2, 
for the Faddeeva function goes to unity in this limit. This ensures, in view of Eq. (|4()jl . 
that lim 4 ^o+ V'( r ;^) = ip{r;t = 0). Equations (J52*Jl and (JBTJj) also allow us to determine the 
long-term behavior of the wave packet in the interval [0, a]. Using the asymptotic expansion 
of the Faddeeva function 0] and Eq. (J57j) . the wave packet for large t is seen to consist of 
bound-state components evolving according to the exp{— iE n t} factor in Eq. (j^S| as well 



as a decaying component that evolves in time as t 3//2 . This speci fic p ower-law 
a decaying state — in the long-time limit — is in fact well known 
0,0] for an alternative view). 



Dehavior of 
see Refs. 



III. NUMERICAL STUDIES 

In this section, we apply the techniques developed so far to a wave packet in the model 
potential of Eq. (JTJ) [V = 5 and a = 10]. For the basis {yj(r) : j = 1, . . . ,N}, a finite-element 



basis set 



BO, E 



3ased on fifth-order Hermite interpolating polynomials was used |2j 



B, 28, 29,0, 



A mesh of evenly-spaced nodes was defined radially, with three finite elements 
centered at the ith node, r^. All three vanish in [0, 7*j_i] and [rj + i,a]. Their behavior at 
separates them into three types. The zero-type polynomials have a finite function value at 
the ith node, but zero first and second derivatives. The one-type have a finite first derivative, 
but zero function value and zero second derivative. The two-type have a finite curvature, 
but vanishing zero and first derivatives. In order to enforce vanishing at the origin, Eq. (JEJ, 
only the one- and two- types are considered at the first node (at r = 0). Thus, the dimension 
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of the finite-element basis set (N) is equal to one less than three times the number of nodes. 

solve the generalized eigenvalue problem of Eq. ()20|) to obtain 2N 
Siegert pseudostates with 2N distinct eigenvalues x n . Appendix |X] proves that certain 
subsets of iV Siegert pseudostates achieve completeness in C^. However, which N states 
must be selected is not clear. By looking at the eigenvalues of the matrix M G <C NxN 
[Eq. ()44|)]. one can determine whether or not a selected subset of N vectors is complete. If 
the vectors form a basis, the M matrix will have N (not necessarily distinct) eigenvalues 
which differ from zero. For the potential under consideration, the number of bound states 
equals the number of antibound states. (In the Introduction, we mentioned that there are 
10 physical bound and nine physical antibound states. The 10th antibound state found in 
the numerical calculation cannot be converged and does not correspond to an eigenstate 
of the Hamiltonian.) All combinations of bound or antibound with incoming or outgoing 
Siegert pseudostates supply iV distinct vectors, giving four plausible choices for a complete, 
useful basis set. The eigenvalues of M were calculated for each case, and in so far as zero 
was never one of them, all combinations proved to span C^. More specifically, more than 
95 % of all eigenvalues of M equal unity to machine precision. The rest have a magnitude 
ranging from 0.93 to 230. 

To assess the ability of the Siegert pseudostates to represent a wave packet at t = 0, it 
is first necessary to investigate the limiting accuracy of the underlying finite-element basis 
set. An initial wave packet is given the form of a Gaussian multiplied by a plane wave: 



The Gaussian is centered in the middle of the well at r = 5, and £ = 0.5 is chosen to ensure 
that contribution outside the interval [0, a] is negligible. 

We represent the wave packet as a superposition of the finite elements, 



ip(r;t = 0) = exp [ 




+ ik (r - r )] . 



(70) 



2£ 2 



N 




(71) 



where the expansion coefficients are chosen to minimize 




(72) 
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i.e., 



ft; 



^(S- 1 ) / yi (r)V>(r;t = 0)dr. 

7 = 1 



(73) 



On the basis of xt, we determine how accurately our finite set of piecewise-defined polyno- 
mials is able to approximate a specific element of the infinite-dimensional Hilbert space. 

Let us next consider the completeness relation expressed by Eq. (|51|). The set of N 
selected Siegert pseudostates must be capable of representing the unique coefficient vector 
a e C N , 

N 

« = ^7n c n, (74) 
n=l 

where, ideally, ct = a. According to Eq. (|51|). 



N 



7* 



y (m- 1 ) cis a . 

/ -i \ ) ran n 



(75) 



n=l 



The quality of the Siegert pseudovector expansion is measured by 



xl 



N 

E 



X, 1 2 



(76) 



Now with the expansion coefficients themselves superpositions of the Siegert pseudostates, 
we calculate 



xl 



N 



ij)(r;t = 0) - ^7n<Pn(r) 



n=l 



dr 



(77) 



to test the ability of the chosen subset of N Siegert pseudostates to accurately reproduce 
the wave packet at t — 0. 

We now focus on Eq. (|40jh which allows us to expand the wave packet in terms of all 2iV 
Siegert pseudostates: 

N 



2 A 



-Ec 



n=l 



Cn — ~ c ,, , 



(78) 
(79) 



xl 



2N 



ip(r;t = 0) - yCnVni 



n=l 



dr . 



(80) 



All five x 2 's were calculated as a function of N, the number of finite-element basis func- 
tions, with k = 15 [see Eq. (|7Dj)]. The values are displayed in Table |U Only x^'s an d x|'s 
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TABLE I: The x 2 's described in Eqs. (1721)- pi . as a function of finite-element basis size, N, 
calculated for the case where k$ = 15 in Eq. (|70l) . xi demonstrates convergence to accurate wave- 
packet reproduction as the number of finite-element basis functions is increased, xi an d x\ show 
the ability of the expansions of Eqs. (|52jl and (|62[). respectively, to match the unique coefficients 
of Eq. 1)73(1 . xi an d X§ confirm that these expansions are successful in reproducing the initial wave 
packet at the level limited by the underlying finite elements. Notation x[y] stands for x x W y . 



N 


xi 


xl 


xl 


xl 


xl 


20 


0.89 


2.7[-19] 


0.89 


9.6[-20] 


0.89 


80 


1.5[-3] 


1.4[-17] 


1.5[-3] 


2.2[-19] 


1.5[-3] 


200 


2.0[-7] 


4.5[-20] 


2.0[-7] 


5.3[-21] 


2.0[-7] 


380 


3.0[-ll] 


7.8[-20] 


3.0[-ll] 


9.8[-21] 


3.0[-ll] 


620 


5.5[-14] 


4.5[-18] 


5.5[-14] 


9.9[-19] 


5.5[-14] 



using bound and outgoing Siegert pseudostates are tabulated, since these states provided 
the qualitatively best results for time evolution based on Eq. ([52)1 . Other combinations of 
N Siegert pseudostates will not be discussed further. 

One observes in Table H] that as the number of basis functions being used is increased, 
Xi approaches convergence. For all values of N, the high accuracy of x\ confirms that the 
Siegert pseudovectors indeed form a basis for . Because of the excellent precision achieved 
at the vector level, we see that xli which measures the representability of the initial wave 
packet in terms of the spatial representation of the N selected Siegert pseudostates, is only 
limited by the degree to which the finite elements are complete. The quantities xl an d x\ 
in Table |J demonstrate that similar statements hold with regard to accuracy if one utilizes 
the overcomplete set of 2N Siegert pseudostates. 

Now consideration is given to the propagation of a wave packet in time. To that end, 
we require a benchmark of comparison. Positive-energy eigensolutions in the interval [0, a] 
were found for the potential of study, Eq. (JTJ), by matching at the discontinuity to energy- 
normalized continuum solutions of the form sin (At + <5). These were used, together 
with the bound-state solutions, to expand the Gaussian wave packet of Eq. f)70j) . We refer 
to this expansion as the "analytical solution," even though it should be pointed out that 
we carried out the integration over the energy continuum numerically. The validity of the 
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FIG. 3: Time evolution of the wave packet described at t = by Eq. (J7UJ); ko = 15. In the r 
interval between and 10, the potential is constant [Vb = 5 and a = 10 in Eq. Iflj)]. 

analytical solution was established by comparison with the closed-form expression in Eq. 
()70|) at t — 0. Numerical convergence of the energy integration was tested carefully and 
ensured. Standard exponential time dependence was introduced for t > 0. 

Once it is determined that the analytical solution gives a wave packet representation 
that is exact to machine precision (at t — 0), we can compare wave packets represented 
by Siegert pseudostates with confidence that discrepancy is due to error on the part of the 
Siegert pseudostates. We can do this using either Eq. (jo^j) or Eq. (jo^j) with Eq. (JBTJj) to 
arrive at expansion coefficients. For both forms we consider two cases: one where the initial 
wave packet has an average energy sufficient to pass the potential step at r = a with little 
reflection (ko = 15), and one where ko = 5, ensuring significant physical reflection. The 
wave packets in the two cases are illustrated in Figs. 01 and |U 

In order to gauge the accuracy of Eq. [i.e., exponential time evolution using a basis 
of iV Siegert pseudostates], we modified xi, Eq- ([77|h to include time dependence. We refer 
to this straightforward extension as xl(t)- An analogous modification of Eq. ()80|) — xi(t)~ 
was introduced to test the nonexponential time evolution expressed by Eqs. (JH2I) and (|o^|). 
Tables |H] and IIHI were calculated with iV = 620. They compare time evolution using Eq. 
(|52~|) with time evolution given by Eqs. (Jfi^ and 4lJ . To account for decaying amplitude 
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FIG. 4: The same parameters as in Fig. |21 but with ko = 5. Notice that there are significant 
physical reflections as soon as the wave packet hits the discontinuity at r = a. 

TABLE II: The time-dependent x 2 ' s discussed in Sec. II I II for the case of the fast-moving wave 
packet (ko = 15; see Fig. x§(^) measures the quality of exponential time evolution, the wave 
packet being expanded in a basis of bound and outgoing Siegert pseudostates [Eq. (|52|)]. x|(t) refers 
to nonexponential time evolution using all 2N Siegert pseudostates [Eqs. (|62[). (|69|) ]. || ?/>(£) || 2 is 
defined in Eq. (|81|), The information in this table confirms that both expansions are accurate in 
describing the motion of a higher-energy particle. 



; Mj xlgv ll gg) ll 2 xl(t)/ ii gffl 

6.2 x 10~ 14 6.2 x 10~ 14 

0.125 1.4 x 1(T 13 1.4 x 1(T 13 

0.25 2.2 x 10~ 13 2.2 x 10~ 13 

0.375 1.5 x 10~ 12 1.5 x 10~ 12 

0.5 1.9 x 10~ 9 1.9 x 10~ 9 



within the region [0,a], we divide X^{t) and xf(t) by 

\\m \?= [ |^M)| 2 dr. (81) 



o 



One sees in Table |H] that for the case where k Q = 15, Eq. (|52|) gives accurate time 
evolution, while Table ITTT1 shows a significant loss of precision for the slower wave packet, 
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TABLE III: The case of the slow-moving wave packet (ko = 5; see Fig. 12). xi(<) tests Eq. (EH> ; 
x|(i) tests Eqs. ()62j) and (|69[). || ip(t) || 2 is defined in Eq. (|81|). For both forms of expansion 
coefficients [Eq. (|52~)l vs. Eq. (|69|)]. the accuracy at t = is only limited by the finite-element 
basis set. For t > 0, xi(^) is many orders of magnitude smaller than x|(t), showing that, for 
this case, accurate time evolution is given by the expansion of Eq. 1)621) using the nonexponential 
time-evolution coefficients in Eq. (|69|) . 



; [a.u.] xi(t)i ii m ii 2 xigv ii m 

2.9 x 10~ 16 2.9 x 10~ 16 
0.5 6.2 x 10~ 6 5.3 x 10~ 16 

1 5.5 x 10" 5 2.2 x 10~ 15 
1.5 4.4 x 10~ 4 2.4 x 10" 14 

2 1.2 x 10~ 3 2.3 x 10~ 13 



ko = 5, at all t > 0. This is because the assumption of pure exponential time evolution 
in Eq. (jo^j) is not justified. In fact, the more complicated time evolution suggested by Eq. 
(|69|) provides superior agreement for both k = 15 and ko = 5. It is important to notice 
that the expansion coefficients associated with all Siegert pseudostates evolve in time in a 
nonexponential fashion. [Equation ()52)) fails also — for k^ = 5 — if Vq is chosen to be so small 
that there are no bound states in the well.] Nevertheless, at least when k$ = 15 (and bound 
and outgoing Siegert pseudostates are selected), Eq. (jo^jl is satisfactory. In this case, we 
found numerically that the vector defined by its elements c^Scx = (ip n \if)), n = 1,...,N, 
is, to machine precision, an eigenvector of the matrix M e c JVxAr with eigenvalue one. As 
a consequence, the coefficients 7 n in Eqs. (|75j) and (jTTjl equal ((p n \ip). This means that 
Eq. (|52j) reduces to Eq. (J3J). We form the conjecture that this is the reason why the 
approach of Refs. |3, Q , based on Eq. (JH), has been successful. [For k = 5, we found 
the performance of Eq. (jlj) to be even worse than that of Eq. (|52|).] The most general and 
accurate treatment, however, is only accomplished on the basis of the time evolution derived 
from the out going- wave Green's function, Eqs. (jr]2^) and 
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IV. CONCLUSION 



In this paper,we reviewed the numerical technique developed by Tolstikhin, Ostrovsky, 
and Nakamura ^, |(| for calculating Siegert pseudostates. The approach is applicable to 
general short-range potentials; it is not restricted to the simple model potential we made 
use of in our numerical investigation. We incorporated a straightforward extension of the 
method to nonorthogonal basis sets and explored in some detail the mathematical properties 
of the Siegert pseudovectors, in particular questions of completeness. 

Our main goal in this paper has been the development of a rigorous formulation of the 
time evolution of a wave packet expanded in terms of Siegert pseudostates. Our study 
demonstrates that assuming pure exponential time evolution can be a poor approximation. 
The technique based on Eqs. (J52j) and ([69)1 . which implies nonexponential time evolution for 
the entire set of 2N Siegert pseudostates, can reproduce very well the analytical reference 
wave packet moving in the model potential. There are no artificial reflections at the boundary 
of the grid (r = a). Physical reflections are correctly reproduced. 

We have shown that it is possible to find a subset consisting of iV Siegert pseudovectors 
that represents a basis of C^. Siegert pseudostates were used in Ref. |3| to carry out 
channel expansions in the context of multichannel quantum defect theory, using bound 
and outgoing Siegert pseudostates only. This paper confirms that such a choice provides 
a complete basis, in principle. In the time domain, this means that a wave packet can be 
represented for any t > as a superposition of these N basis states. Their time evolution 
is, however, nonexponential. The time-dependent expansion coefficients may be found by 
applying the completeness relation implied by Eq. (|51|) to the expansion in Eq. (|62|) . 

While wave-packet propagation using Eqs. (jo^j) and (JoTJj) is numerically stable and ac- 
curate, the calculation of all 2N eigenvalues and eigenvectors of the generalized eigenvalue 
problem in Eq. ()20|) is necessary in principle. This requirement poses a challenge for itera- 
tive sparse-matrix techniques. Further developments will be needed in order to turn Siegert 
pseudostates into a useful tool for large-scale calculations on atoms and molecules. Siegert 
pseudostates have already offered a glimpse at their extraordinary potential. 
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APPENDIX A: EXISTENCE OF A SUBSET OF THE 2N SIEGERT PSEU- 
DOVECTORS FORMING A BASIS OF 

Let Nq denote the number of eigenvectors c n whose eigenvalue x n satisfies Im(x n ) > 
[Re(k n ) > 0]. There is an identical number of eigenvectors c n with Im(x n ) < [Re(A; n ) < 0]. 
The remaining 2N — 2Nq eigenvectors are characterized by Im(x n ) = [Re(fc n ) = 0]. The 
latter correspond to the bound and the antibound states. We order the 2N eigenvectors as 
follows: 

Ci , . . ■ , C Nq j Cjy c+ i , . . . , C2N-N C j C 2 N-N C +U ■ ■ ■ ■> c 2jV • (Al) 
lm(j< n )>0 lm(>c n )=0 lm(*»)<0 

Additionally, the vectors c 2 n-n c +1i ■ • ■ > C 2N are ordered in such a way that 

c 2N -N c +n = c* , n = 1, . . . , N c (A2) 

[note the remarks following Eq. We now describe a constructive approach for ex- 

pressing the last Nq Siegert pseudovectors in Eq. (jAlj) [Re(k n ) < 0] in terms of the vectors 

Ci, . . . , C 2 N-N C - 

Equations and pUj) can be multiplied from the right by Sc m . Let, in particular, 
1 < m < Nc- Then we can exploit that c^Scm > and solve for c^: 

2N ~ Nc x* r T <?r 



- E («> 

1=1 

N C 

E 



! C m Sc m 

V " xIclSc 



n=l(n^m) 



2^V C To r 



„=l(„^ m ) c ™*c m 
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According to these equations, c 2 n = c*n c can be represented in terms of c±, . . . , c 2 n-i, i-e., 

rank([ci, . . . , c 2N -i]) = rank([ci, . . . , c 2N \) = N . (A5) 



Let us assume now that, using Eqs. (jA3|) and (jA4|) . we have been able to show for a fixed 
m, 1 < in < Nq, that 

2N~N C N c ~m 

c* Nc _ n+1 = a fc° c * + E Pk l)c k ,n = l,...,m (A6) 

k=l k=l 

for suitably chosen a^\j3^ G C. This is clearly true for m = 1 [Eq. (|A5|) ]. It follows from 
Eq. ()A6|) that rank(C) = N, where C = [ci, . . . , C2jv-m]- Also note that the coefficients 
a h ] ifik^ cannot be unique (since 2N — m > N), unless m = Nq and Nc = N. In fact, 
there are an infinite number of solutions. This is easy to see: Any solution e <^N-m 
of CzW = c* Nc _ n+ i can be written as the sum of a particular solution of this linear system 
(which exists, by assumption) and a solution y e t^ 2N - m Q f the homogeneous system, Cy = 
0. The kernel of the matrix C is nontrivial, for its dimension is dim(C 2iY ~ m ) — rank(C) = 
N — m. Hence, there are infinitely many vectors y satisfying the homogeneous system. 
We must make the step from mtom + 1. It is sufficient to demonstrate that 

2N-N C Nc-m-1 

c*N C - m = E «t m+1) ^+ E ^ m+1) <- ( A? ) 

k=l k=l 

Equation (jA3|) . applied to c* N , can be written as 



2N-N C Nc-m-1 

-Nc~m 



E (•••K+ E (■■■)< ( A§ ) 



n=l n=l 



\ -» X V c -m C, ^Cjy c _ m , 



n=V c -m+l X ™ C N c -m^ C Nc-m 



Combined with Eq. (jA6|) . this yields 



2N-N C Nc-m-1 

*N c -m = E (■ ■ -) C ™ + E (" • 0< 

n=l n=l 



n= N c -m+l X ™ C N c -m^ C ^C-m 
X @Nc-m C N c -m • (A9) 
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The exact form of the coefficients symbolized by (...) is inessential. If we use Eq. (jA4j) in 
place of Eq. (jA3|) , we find 

2N-N C N c -~m-l 

c* Nc _ m = Yl (•••)c„+ J2 

n=l n=l 

N c + C 

\ ^ C n OCN c _ m 



n= N c -m+l C N c -m^ C N C -m 
X @Nc-m + C *N c -m ■ (A10) 



Therefore, as long as either 



N C * f c 



n= N c -m+l n '-N c -m u -"-Nc-m 

or 

N c t c 

E C n OC N c -m n(N c -n+l) 

n=N c -m+l ^Nc-m^^Nc-m 



differ from one, Eq. (|A7J) is satisfied. Otherwise, since the coefficients (3^ are not unique, 

'N c -m 



it is possible in general to choose the /3^ < i r ^ +1 ' > such that either Eq. (jA9|) or ()A10|) can be 



solved for c* Nc _ m . 

The coefficients /3^ c L ri ^ +1 ' ) will be unique, however, if and only if 

rankQci, . . . , c 2 jv_ m ]) = rank([ci, . . . , c 2 jv-m-i]) + 1 , (All) 

which would also be consistent with the disappearance from Eqs. (jA9|) and (jA10|) of the 
term involving c* Nc _ m . Thus, in exact arithmetic, it is conceivable that the construction 
step from mtom + 1 fails, i.e., we cannot prove that all incoming Siegert pseudovectors can 
be eliminated. (Nevertheless, there is no reason to anticipate this to cause any difficulties 
in numerical calculations.) If really necessary, then the vector c* Nc _ m must be kept as an 
essential basis vector, and the construction procedure outlined above may be continued for 
the remaining vectors c 1; . . . , c 2 N-m-i- 

Not only the vectors c* [Re(/c n ) < 0] can be eliminated (except for the possible — though 
unlikely — occurrence of essential vectors), but equations analogous to Eqs. (|A3J) and (|A4|) 
exist also for the bound and the antibound eigenvectors. Note that for these, c^Scm differs 
from 0. It is therefore clear that linearly independent basis vectors c n , n = 1, . . ., N, 
can be found among the vectors c±, . . . , C2N- The completeness relation satisfied by these A" 
vectors is given by Eq. (pTTj) . 
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APPENDIX B: TIME EVOLUTION AND THE FADDEEVA FUNCTION 



For evaluating the integral in Eq. (|66|). it is useful to perform the substitution 

s 2 = ikH/2 . (Bl) 

Thus, 



"oo „-ik 2 t/2 r+y/it/2oc -s 2 



dk = I _ 6 ds . (B2) 



, k - k n J-^/u/2oo s - ^'\tj1k n 
The right-hand side of this equation is reminiscent of the Faddeeva function defined in Eq. 
(IfjHjl . However, the integration path must be rotated back to the real axis. In order to do 
this, we must carefully distinguish between poles due to bound, antibound, outgoing, or 
incoming Siegert pseudostates. The pole of the integrand on the right-hand side of Eq. ()B2jl 
lies in the upper s plane if, for t > (v/i/2 > 0), 



I = eW1 vl < B3) 

for bound or outgoing Siegert pseudostates, and 

| = ^"V! (B4) 

for antibound or incoming Siegert pseudostates. (Note that the Faddeeva function is defined 
in the upper half-plane.) 

Let us first consider a bound state. In this case, the pole is in the second quadrant of 
the complex s plane. The integration path, illustrated in Fig. together with the location 
of the bound-state pole, runs along the diagonal from the third to the first quadrant. Since 
the contour may be closed as indicated in the figure, and since there are no poles inside the 
resulting loop, we find 

/it/ioo e -s 2 

, — ds = 

-^it/200 S - yj\t/2k n 




(B5) 



for a bound state. 
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Im(s) 




Re(s) 



lm(s) 





lm(s) 




FIG. 5: The integration contour used to evaluate Eq. (|B2|) depends on the nature of the Siegert 
pseudostate considered. Arrows indicate the direction of the integration path in the complex s 
plane. The circle symbolizes the location of the pole of the integrand on the right-hand side of Eq. 
I)B2|1 . a) Bound state, b) Outgoing Siegert pseudostate. c) Antibound state, d) Incoming Siegert 
pseudostate. 

The location of the pole in the case of an outgoing Siegert pseudostate is shown in Fig. 
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03d. This time the pole lies inside the integration loop, so that 

-ds = -2me-' lk " t/2 (B6) 



'it/2oo Q - s 2 



-Jit/200 s - y/it/2k r . 



■V" Ws* 



n 



for an outgoing Siegert pseudostate. 

The pole's location for an antibound and an incoming Siegert pseudostate, respectively, 
can be seen in Figs. and (The integration path now runs along the diagonal from 
the first to the third quadrant.) In neither case does the pole interfere with the rotation 
of the integration path to the real s axis. Therefore, for an antibound or incoming Siegert 
pseudostate the result reads 

y/it/2oo g— s 2 poo Q~ s2 



ds = — ■ ; — — ds 



I^ioo s - ^Jttj2k n J-oo s + eP^y/tftk 



-e^/y ^ I . (B7) 
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